# yI1-xI0-EMWproposal.R
# Philips, Andrew Q. "How to Avoid Incorrect Inference (While Gaining Correct Ones) in Dynamic Models". Forthcoming at Political Science Research and Methods. 
# andrew.philips@colorado.edu
# 2/9/21
#
# Figures produced in this R script:
#   --SI, Figure 32: "yi1-xi0-LRlxYES-EMWproposal.pdf" 
# -----------------------------------

# -----------------------------------
# EXPERIMENT 3: Y ~ I(1), X ~ I(0)
# Type I error
# ONLY CALCULATING THE LRE IF L.X (ARDL/ECM) AND X (LDV) ARE FIRST STAT SIG (Section 2)
# -----------------------------------
rm(list = ls())
setwd("~/Dropbox/My_Folder/Papers-Projects/Wlezien-Enns-PSRM/Final-Feb2021/scripts/")
# -----------------------------------
library(dynamac) 
library(nlWaldTest)
library(ggplot2)
library(grid)
library(gridExtra)
library(dplyr)

set.seed(102938475)
# -------------------------------------------------------------
# EXPERIMENT 3: Y ~ I(1), X ~ I(0)
sims = 2000 # no. sims
ar.param <- seq(0,0.95, by = 0.05) # autoregression in X
N <- c(50, 250, 1000) # number of obs
e.var <- 1  #c(1, 5) # error variance in Y

container.ardl <- matrix(NA, nrow = sims*length(ar.param)*length(N)*length(e.var), ncol = 10) # null matrix to hold output
# cols = sims, ar.param.X, N, Beta.lx, UL.beta.lx, LL.beta.lx, LR.x, LR.x.UL, LR.x.LL, e.var
container.ecm <- matrix(NA, nrow = sims*length(ar.param)*length(N)*length(e.var), ncol = 10) # null matrix to hold output
# cols = sims, ar.param.X, N, Beta.lx, UL.beta.lx, LL.beta.lx, LR.x, LR.x.UL, LR.x.LL, e.var
container.ldv <- matrix(NA, nrow = sims*length(ar.param)*length(N)*length(e.var), ncol = 10) # null matrix to hold output
# cols = sims, ar.param.X, N, Beta.x, UL.beta.x, LL.beta.x, e.var
# -------------------------------------------------------------

row <- 1
prog <- txtProgressBar(min = 0, max = sims, style = 3)
for (s in 1:sims) {
  for (ar.x in ar.param) {
    for (n in N) {
      for (ev in e.var) {
        setTxtProgressBar(prog, value = s)
        # create X
        if (ar.x == 0) {
          x <- rnorm(n + 100)
        } else {
          x <- arima.sim(n = (n + 100), list(ar = c(ar.x)))
        }
        # create Y
        y <- cumsum(rnorm(n + 100, sd = sqrt(ev)))
        # ARDL model:
        res.ardl <- lm(y[101:length(y)] ~ lshift(y, 1)[101:length(y)] + x[101:length(x)] + lshift(x, 1)[101:length(x)])
        # ECM model:
        res.ecm <- lm(dshift(y)[101:length(y)] ~ lshift(y, 1)[101:length(y)] + dshift(x)[101:length(x)] + lshift(x, 1)[101:length(x)])
        # LDV model
        res.ldv <- lm(y[101:length(y)] ~ lshift(y, 1)[101:length(y)] + x[101:length(x)])
        # static model
        res.static <- lm(y[101:length(y)] ~ x[101:length(x)])
        
        # grab QOI (FORALL): ----
        container.ardl[row, 1] <- container.ecm[row, 1] <- container.ldv[row, 1] <- s
        container.ardl[row, 2] <- container.ecm[row, 2] <- container.ldv[row, 2] <- ar.x
        container.ardl[row, 3] <- container.ecm[row, 3] <- container.ldv[row, 3] <- n
        container.ardl[row, 10] <- container.ecm[row, 10] <- container.ldv[row, 10] <- ev
        
        # QOI (ARDL)
        container.ardl[row, 4] <- coef(res.ardl)[4] # beta.lx
        container.ardl[row, 5] <- confint.default(res.ardl, parm = c(4), level = 0.95)[2] # UL
        container.ardl[row, 6] <- confint.default(res.ardl, parm = c(4), level = 0.95)[1] # LL
        # LR effects (ARDL)
        lrm <- nlConfint(res.ardl, "(a[3] + a[4])/(1-a[2])", level = 0.95)
        container.ardl[row, 7] <- lrm[1] # LRM
        container.ardl[row, 8] <- lrm[3] # LRM UL
        container.ardl[row, 9] <- lrm[2] # LRM LL
        
        # QOI (ECM)
        container.ecm[row, 4] <- coef(res.ecm)[4] # beta.lx
        container.ecm[row, 5] <- confint.default(res.ecm, parm = c(4), level = 0.95)[2] # UL
        container.ecm[row, 6] <- confint.default(res.ecm, parm = c(4), level = 0.95)[1] # LL
        # LR effects (ECM)
        lrm <- nlConfint(res.ecm, "(a[4])/(-a[2])", level = 0.95)
        container.ecm[row, 7] <- lrm[1] # LRM
        container.ecm[row, 8] <- lrm[3] # LRM UL
        container.ecm[row, 9] <- lrm[2] # LRM LL
        
        # QOI (LDV)
        container.ldv[row, 4] <- coef(res.ldv)[3]
        container.ldv[row, 5] <- confint.default(res.ldv, parm = c(3), level = 0.95)[2] # UL
        container.ldv[row, 6] <- confint.default(res.ldv, parm = c(3), level = 0.95)[1] # LL
        # LR effects (LDV)
        lrm <- nlConfint(res.ldv, "(a[3])/(1-a[2])", level = 0.95)
        container.ldv[row, 7] <- lrm[1] # LRM
        container.ldv[row, 8] <- lrm[3] # LRM UL
        container.ldv[row, 9] <- lrm[2] # LRM LL
        
        row <- row + 1
      }
    }
  }
}
#save.image("scenario-yi1-xi0-EMWproposal-data.RData")
#load("scenario-yi1-xi0-EMWproposal-data.RData")

# Create rejection rates and label:
container.ardl <- as.data.frame(container.ardl)
# cols = sim, ar.param.X, N, Beta.lx, UL.beta.lx, LL.beta.lx, LR.x, LR.x.UL, LR.x.LL, "e.var"
colnames(container.ardl) <- c("sims", "ar.param.X", "Time", "Beta.lx", "Beta.lx.ul", "Beta.lx.ll", "LR.x", "LR.x.UL", "LR.x.LL", "e.var")
container.ecm <- as.data.frame(container.ecm)
# cols = sim, ar.param.X, N, Beta.lx, UL.beta.lx, LL.beta.lx, LR.x, LR.x.UL, LR.x.LL, "e.var"
colnames(container.ecm) <- c("sims", "ar.param.X", "Time", "Beta.;x", "Beta.lx.ul", "Beta.lx.ll", "LR.x", "LR.x.UL", "LR.x.LL", "e.var")
container.ldv <- as.data.frame(container.ldv)
# cols = sim, ar.param.X, N, Beta.x, UL.beta.x, LL.beta.x, LR.x, LR.x.UL, LR.x.LL, "e.var"
colnames(container.ldv) <- c("sims", "ar.param.X", "Time", "Beta.x", "Beta.x.ul", "Beta.x.ll", "LR.x", "LR.x.UL", "LR.x.LL", "e.var")

# proportion false rejection
container.ardl$reject.beta.lx <- ifelse(container.ardl$Beta.lx.ll > 0 | container.ardl$Beta.lx.ul < 0, 1, 0)
container.ecm$reject.beta.lx <- ifelse(container.ecm$Beta.lx.ll > 0 | container.ecm$Beta.lx.ul < 0, 1, 0)
container.ldv$reject.beta.x <- ifelse(container.ldv$Beta.x.ll > 0 | container.ldv$Beta.x.ul < 0, 1, 0)

# coverage of LRM:
container.ardl$reject.lrm <- ifelse(container.ardl$LR.x.LL > 0 | container.ardl$LR.x.UL < 0, 1, 0)
container.ecm$reject.lrm <- ifelse(container.ecm$LR.x.LL > 0 | container.ecm$LR.x.UL < 0, 1, 0)
container.ldv$reject.lrm <- ifelse(container.ldv$LR.x.LL > 0 | container.ldv$LR.x.UL < 0, 1, 0)

# Coding of rejection of LRM:
# = 1 if we reject LRM AND also reject beta.lx/beta.x, 0 otherwise
container.ardl$reject.lrm.yes.lx <- ifelse(container.ardl$reject.lrm == 1 & container.ardl$reject.beta.lx == 1, 1, 0)
container.ecm$reject.lrm.yes.lx <- ifelse(container.ecm$reject.lrm == 1 & container.ecm$reject.beta.lx == 1, 1, 0)
container.ldv$reject.lrm.yes.x <- ifelse(container.ldv$reject.lrm == 1 & container.ldv$reject.beta.x == 1, 1, 0)

# = 1 if we reject LRM AND NOT reject beta.lx/beta.x, 0 otherwise
container.ardl$reject.lrm.no.lx <- ifelse(container.ardl$reject.lrm == 1 & container.ardl$reject.beta.lx == 0, 1, 0)
container.ecm$reject.lrm.no.lx <- ifelse(container.ecm$reject.lrm == 1 & container.ecm$reject.beta.lx == 0, 1, 0)
container.ldv$reject.lrm.no.x <- ifelse(container.ldv$reject.lrm == 1 & container.ldv$reject.beta.x == 0, 1, 0)

# collapse everything down:
container.collapse.ardl <- container.ardl %>%
  group_by(ar.param.X, Time, e.var) %>%
  summarize(reject.beta.lx = mean(reject.beta.lx),
            reject.lrm = mean(reject.lrm),
            reject.lrm.yes.lx = mean(reject.lrm.yes.lx),
            reject.lrm.no.lx = mean(reject.lrm.no.lx))

container.collapse.ecm <- container.ecm %>%
  group_by(ar.param.X, Time, e.var) %>%
  summarize(reject.beta.lx = mean(reject.beta.lx),
            reject.lrm = mean(reject.lrm),
            reject.lrm.yes.lx = mean(reject.lrm.yes.lx),
            reject.lrm.no.lx = mean(reject.lrm.no.lx))

container.collapse.ldv <- container.ldv %>%
  group_by(ar.param.X, Time, e.var) %>%
  summarize(reject.beta.x = mean(reject.beta.x), 
            reject.lrm = mean(reject.lrm),
            reject.lrm.yes.x = mean(reject.lrm.yes.x),
            reject.lrm.no.x = mean(reject.lrm.no.x))

# separate by T:
container.collapse.ardl.50 <- subset(container.collapse.ardl, Time == 50)
container.collapse.ardl.250 <- subset(container.collapse.ardl, Time == 250)
container.collapse.ardl.1000 <- subset(container.collapse.ardl, Time == 1000)
container.collapse.ecm.50 <- subset(container.collapse.ecm, Time == 50)
container.collapse.ecm.250 <- subset(container.collapse.ecm, Time == 250)
container.collapse.ecm.1000 <- subset(container.collapse.ecm, Time == 1000)
container.collapse.ldv.50 <- subset(container.collapse.ldv, Time == 50)
container.collapse.ldv.250 <- subset(container.collapse.ldv, Time == 250)
container.collapse.ldv.1000 <- subset(container.collapse.ldv, Time == 1000)

# plot: -------------------------------------------
# LR, T = 50
p1 <- ggplot() + 
  geom_line(data = subset(container.collapse.ardl.50, e.var == 1), aes(x = ar.param.X, y = reject.lrm.yes.lx), color = "#a63603", size = 1.7, lty = 'solid') + 
  geom_line(data = subset(container.collapse.ecm.50, e.var == 1), aes(x = ar.param.X, y = reject.lrm.yes.lx), color =  "#e6550d", size = 1.7, lty = 'dotted') + 
  geom_line(data = subset(container.collapse.ldv.50, e.var == 1), aes(x = ar.param.X, y = reject.lrm.yes.x), color = "#fd8d3c", size = 1.7, lty = 'dashed') + 
  scale_x_continuous(breaks = seq(0,.95,.15)) +  scale_y_continuous(limits = c(0, 0.4)) + geom_hline(yintercept = 0.05, color = "black", size = 1) +
  ylab("Proportion Rejected") + xlab("Autoregression in X") + ggtitle("Long-Run, T=50") +
  theme(panel.background = element_rect(fill = NA), panel.grid.major = element_line(colour = "grey88"), panel.ontop = FALSE, # customize axis lines in background
        text = element_text(size=10), # size of axis numbers
        plot.title = element_text(size = 15)) 

# LR, T = 250
p2 <- ggplot() + 
  geom_line(data = subset(container.collapse.ardl.250, e.var == 1), aes(x = ar.param.X, y = reject.lrm.yes.lx), color = "#a63603", size = 1.7, lty = 'solid') + 
  geom_line(data = subset(container.collapse.ecm.250, e.var == 1), aes(x = ar.param.X, y = reject.lrm.yes.lx), color =  "#e6550d", size = 1.7, lty = 'dotted') + 
  geom_line(data = subset(container.collapse.ldv.250, e.var == 1), aes(x = ar.param.X, y = reject.lrm.yes.x), color = "#fd8d3c", size = 1.7, lty = 'dashed') + 
  scale_x_continuous(breaks = seq(0,.95,.15)) +  scale_y_continuous(limits = c(0, 0.4)) + geom_hline(yintercept = 0.05, color = "black", size = 1) +
  ylab("Proportion Rejected") + xlab("Autoregression in X") + ggtitle("Long-Run, T=250") +
  theme(panel.background = element_rect(fill = NA), panel.grid.major = element_line(colour = "grey88"), panel.ontop = FALSE, # customize axis lines in background
        text = element_text(size=10), # size of axis numbers
        plot.title = element_text(size = 15)) 

# LR, T = 1000
p3 <- ggplot() + 
  geom_line(data = subset(container.collapse.ardl.1000, e.var == 1), aes(x = ar.param.X, y = reject.lrm.yes.lx), color = "#a63603", size = 1.7, lty = 'solid') + 
  geom_line(data = subset(container.collapse.ecm.1000, e.var == 1), aes(x = ar.param.X, y = reject.lrm.yes.lx), color =  "#e6550d", size = 1.7, lty = 'dotted') + 
  geom_line(data = subset(container.collapse.ldv.1000, e.var == 1), aes(x = ar.param.X, y = reject.lrm.yes.x), color = "#fd8d3c", size = 1.7, lty = 'dashed') + 
  scale_x_continuous(breaks = seq(0,.95,.15)) +  scale_y_continuous(limits = c(0, 0.4)) + geom_hline(yintercept = 0.05, color = "black", size = 1) +
  ylab("Proportion Rejected") + xlab("Autoregression in X") + ggtitle("Long-Run, T=1000") +
  theme(panel.background = element_rect(fill = NA), panel.grid.major = element_line(colour = "grey88"), panel.ontop = FALSE, # customize axis lines in background
        text = element_text(size=10), # size of axis numbers
        plot.title = element_text(size = 15)) 

pdf("yi1-xi0-LRlxYES-EMWproposal.pdf", width = 1.618*7, height = 7)
grid.arrange(p1, p2, p3, ncol = 3)
dev.off()
